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A procedure is introduced for deriving a coarse-grained dissipative particle dynamics from molec- 
ular dynamics. The rules of the dissipative particle dynamics are derived from the underlying molec- 
ular interactions, and a Langevin equation is obtained that describes the forces experienced by the 
dissipative particles and specifies the associated canonical Gibbs distribution for the system. 
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Hydrodynamic simulations of complex fluids remain a 
major challenge in most cases of interest. Such fluids 
include particulate and colloidal suspensions, polymeric 
liquids, emulsions and other self-assembling amphiphilic 
fluids, and fluids where Brownian motion is important. 
For such fluids it is often necessary to base the model- 
ing on a microscopic picture of the system, thus working 
from the bottom upwards. Over the last decade several 
such 'bottom up' strategies have been introduced. Hy- 
drodynamic lattice gases ||^, which model the fluid as a 
discrete set of particles, represent a computationally ef- 
flcicnt discretization of the more conventional molecular 
dynamics (MD) [||. 

A recent contribution to the family of bottom-up 
approaches is the dissipative particle dynamics (DPD) 
method introduced by Koelman and Hoogerbrugge in 
1992 HJ. Applications of the technique include colloidal 
suspensions Q , polymer solutions and binary immisci- 
ble fluids ^ . For speciflc applications where comparison 
is possible, this model is orders of magnitude faster than 
. MD §. 

The basic components of DPD are particles that are 
thought to represent mesoscopic elements of the under- 
lying molecular fluid. These dissipative particles then 
evolve just as MD particles but with different inter- 
particle forces: Since the DPD particles are pictured as 
having internal degrees of freedom, the forces between 
them have both a fluctuating and a dissipative com- 
ponent in addition to the conservative forces that are 
present already at the MD level. Nevertheless, momen- 
' tum conservation along with mass conservation produce 
hydrodynamic behavior at the macroscopic level. 

Dissipative particle dynamics has been demonstrated 
to connect correctly to the macroscopic continuum the- 
ory; that is, for a one-component DPD fluid, it is possible 
to derive the Navier-Stokes equations and to compute the 
viscosity in the large scale limit [^,^. However, thus far 
no attempt has been made to link DPD to the under- 
lying microscopic dynamics. This is the purpose of the 
present letter. We deflne the dissipative particles (DP) 



by appropriate weight functions that sample a portion 
of the underlying conservative MD particles, and we de- 
rive the forces between the DP's from the hydrodynamic 
description of the MD system. 




FIG. 1. Multiscale modeling: The dissipative particles 
are defined as cells in the Voronoi lattice, moving with veloc- 
ity Ufe. There are four relevant length scales: The scale of the 
large, gray colloid particles, the two scales of the dissipative 
particles in between and away from the colloids and finally 
the scale of the MD particles, which are shown as the little 
dots that form the boundaries between the DP's. 

The present development has two main virtues, one 
fundamental and one practical. From a fundamental 
point of view our work gives a microscopic foundation 
to DPD and thus provides a quantitative meaning to the 
term 'mesoscopic'. On the practical side this foundation 
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may be used to deal with physical systems where the 
modeling is challenged by the simultaneous presence of 
several different length scales. While conventional DP's 
are spheres of fixed size and mass, the current DP's are 
defined as cells on a Voronoi lattice with variable sizes. 
This provides us with the freedom to define particle sizes 
according to the local resolution requirements-a particle 
analog to adaptive meshes in finite-element simulations 
I pof . The concept is illustrated by the simulation of a 
colloidal suspension, which is shown in Fig. Here the 
computational effort is adapted to meet the local need 
for detail of description: it is larger in narrow regions be- 
tween the particles than in the bulk. Previous DPD sim- 
ulations have had difficulty with dense colloidal suspen- 
sions precisely because the technique is unable to handle 
multiple lengthscale phenomena . Other complex sys- 
tems where modeling and simulation frequently involve 
several simultaneous length scales include polymeric and 
amphiphilic fluids, particularly in porous media and re- 
stricted geometries pl|] . 

The basic ingredient in our derivation of DPD is an ap- 
propriate coarse-graining scheme. The dissipative parti- 
cles are defined as clusters of MD particles in such a way 
that the MD particles are all represented by the dissipa- 
tive particles. A general way to achieve this is via the 
sampling function /fc(x) = s(x — Vk)/'}2,i — r/). Here 
the positions and r; define the DP centers, which ini- 
tially may be distributed arbitrarily in space, x is an 
arbitrary position and s(x) is some localized function, 
which we choose as a Gaussian s(x) — exp {—x^ / a?'); the 
distance a sets the scale of the sampling function. The 
mass Mfc, momentum and internal energy Ek of the 
fcth DP are then defined as 



Mk = y^/fc(x^)m 

i 

i 
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ter some algebra the microscopic equations of motion 
then take the form 



(1) 



where x^ and v.; are the position and velocity of the 
i'th MD particle, all assumed to have identical masses 
m, V{rij) is the MD interparticle potential, and is 
the velocity of the fc'th DP. The kinematic condition 
ffc = Ufc completes the definition of the DPD. The nor- 
malization property X]fc/fe(^) = 1 implies directly that 
J2k^^k = and X]fc^'4Ufc = Z]i™Vi, so that if 

mass, momentum and energy are conserved at the MD 
level, they are also conserved at the DP level. 

In order to obtain the equations of motion for the 
DPD we now take the time derivatives of Eqs. (|l|). The 
Gaussian form of s makes it possible to write the time- 
derivative /fe(xi) = /fc;(xj)(v^ ■ Vki + x^ • Ufej) where the 
function is defined as /fci(x) = (2/a2)/fc(x)/;(x). Af- 



dt 



^ = ^ Mki = Ai(x.)m(v^ • ra + x^ ■ V^i) 



I 



dt 



-JT = E ^ Hr +2^^l{^^) J* - n,r — • vki (2) 



where we have defined the general momentum-flux tensor 
II- = mv^v- -I- (1/2) Yl,j Axy , where F^j is the force 
between MD particles i and j, and the microscopic en- 
ergy flux vector = e.vH(l/4) E.^j '(v'+v^OAxy . 
We have also used the definitions v- = v.^ — (Ufc + U;)/2, 
X- = Xt - (rfc + r;)/2, Uki = Uk - U;, rfe; = rfe - r; 
and mg is the external force on an MD particle. In the 
mass equation above the x^ ■ Ufc; term may be shown to 
be negligible upon averaging, as it samples the difference 
in mass density rather than the average of this quantity 
across the region where fki ^ |I2|. For that reason it 
will be omitted, and we have already omitted the corre- 
sponding terms in the momentum and energy equations. 



A. 
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FIG. 2. The overlap region between two Voronoi cells 
is shown in grey. The sampling function /^(r) is shown in 
the upper graph and the overlap function in the lower graph. 
The width of the overlap region is a^/\rk — r;| and its length 
is denoted by Iki- 

All the interaction terms in the above trans- 
port equations are weighted by the overlap func- 
tion fki{x). If only two DP's, k and I say, 
are present it may be shown that /fc;(x) = 
(f /(2a2)) cosh-2 ((x - (rfc + ri)/2) ■ (r^ - vi)/a^). This 
function becomes exponentially small away from the di- 
viding line that is equally far from r^. and rj, as is il- 
lustrated in Fig. 1^. The set of all such dividing lines 
defines a Voronoi lattice. In Fig. ^ fictitious MD parti- 
cles are plotted where /fci(x) > 0.2a^. This happens in 
the neighborhood of the dividing lines. When additional 
DP's are present their contribution to the cosh"^ result 
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may be shown to be negligible in the vicinity of the divid- 
ing line, except at the corners, where dividing lines meet. 
In the end the DPD equations of motion will turn out 
to be independent of a, and only the length l^i shown 
in Fig. 1^ will remain. At this point it suffices to con- 
struct the Voronoi lattice itself, and there is no need to 
evaluate the overlap functions. Standardized algorithms 
and software for the construction of Voronoi lattices exist 

@- 

Note that since the right hand side terms of the mass 
and momentum part of Eqs. are all odd under the 
exchange I ^ k the DP's interact in a pairwise and ex- 
plicitly conservative fashion. The same is true for the 
energy equation if it is rewritten in terms of the total, 
rather than the internal energy. 

Splitting Eqs. (||) into fluctuating and average parts 
gives 



dMk 
dt 

dPfc 
dt 



^ fki (xi)m(v-) -rki+Y^ Mki 
li I 



^ = ^/..(x.)((J'.)-(n:).^).r 



■t ki ■ — H Qkl ■ 

I 
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where Mki is the fluctuating part of the mass flux, Fki is 
the fluctuating part of the momentum flux fki (xi)n' • 
rki + Mki(^k + Vi)/2, and qki the fluctuating part of the 
energy flux J2^ /fci(x,)J^ • rki + (l/2)Mfc,(Ufc;/2)2. Note 
that we have absorbed the contributions from the mass 
variations in F^; and qm- The thermal averages, (...), 
are computed by means of an ensemble of systems with 
common instantaneous values of the mesoscopic variables 
{rfc, Mk, Ufc, Ek}. This means that only the time deriva- 
tives of this set have a fluctuating part. 

It is necessary to introduce some average description 
of the MD system. For this purpose we assume the scale 
separation a << Irj. — r;|, for all k and I, and that the 
width of the overlap region /rki is larger than the mean 
free path of the MD particles. For simplicity we choose 
the momentum flux tensor of a simple Newtonian fluid 
which has the form (H^;) — IP — ?7(Vv -I- (Vv)-^) where 
77 is the dynamic viscosity and P the pressure of the MD 
fluid, ^ denotes the transpose and I is the identity tensor 
. We shall make the approximation that the average 
molecular particle velocity (v) interpolates linearly be- 
tween the DP's. 

It follows that the average mass current (v') =0 and 
that the velocity gradients in the momentum-flux tensor 
take the form Vv -I- (Vv)^ = l/rki{ekilJki + VkiGki), 
where e^; = (r^ — ri)/\rk — ri\. Since (v') « we may 
choose the heat flux according to Fouriers law (J') = 



AVT, where T is the temperature and A is the thermal 
conductivity. In other words, in the frame of reference of 
the overlap region the energy flux is simply the heat flux 
since work terms proportional to the velocity vanish [ p^ . 

With this input we get Mk = J2i ^ki and 
—— = Mkg - > Iki -rreki + — (Uk 

dt ^-^ \ 2 r,., 



Tkl 



'kl 



(Ufe; • eki)eki) 



dEk 
dt 



+ E^ 

I 

-E^ 



kl 



Pk +Pl 

Ik I — t: — ^kl 



— {Uki + (Ufci • eki)eki)] ■ 
rki J 2 



E' 



Tki 
h.kX — 



U 



kl 



Tkl 



' kl 



Qkl 



where we have assumed that the pressure p and temper- 
ature T, as well as the average velocity, interpolates lin- 
early between DP centers, pki = Pk —pi and Tki = Tk—Ti. 
The pressure will eventually follow from an equation of 
state of the form pk ~ p{Ek,Vk) where Vk is the vol- 
ume of DP k. The pressure and temperature must be 
obtained via a thermodynamic description, i.e. an equa- 
tion of state that relates pressure and temperature to 
energy Ek and volume Vk ■ In the special case of an ideal 
gas (see Ref. for a more general treatment), these 
relations simplify to dpk — Pk(dEk/Ek — dVk/Vk) and 
dTk = dEkl{NkB). 

The average rate of change of Mk vanishes. This allows 
us to neglect mass variations altogether in the DPD equa- 
tions, since the effect of mass fluctuations may then be 
absorbed in F^; and qki- Had there been a coupling, say, 
between the averaged momentum and mass values the 
mass would have had to be updated as well. With noth- 
ing but fluctuations in Mk the only change introduced 
by the Mk =:const. approximation is the loss of the drift 
in the M^'s around their constant average, caused by the 
fluctuations. Nothing is neglected in the instantaneous 
changes of momentum and energy. 

In general, force fluctuations will cause mass fluctua- 
tions, which in turn will couple back to cause momentum 
fluctuations. The time scale over which this will happen 
is tri = ffe;/'7, where 77 is the dynamic viscosity of the 
MD system. This is the time it takes for a velocity per- 
turbation to decay over a distance rki- We shall need 
to make the assumption that the fluctuating forces are 
Markovian, and it is clear that this assumption may only 
be valid on time scales larger than t,,. Since the time 
scale of a hydrodynamic perturbation of size say, is 
also given as P /rj this restriction implies the scale sepa- 
ration requirement r^,; « consistent with the scale 
rki being mesoscopic. 

In 'conventional' DPD P,^]]!^] the forces are pairwise 
and act parallel to Bki- They have a conservative part 
that depends only on rki and a dissipative part propor- 
tional to (Ufei ■ eki)eki- Here the same terms are present. 
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The conservative force is seen to arise from the pressure 
and the dissipative part from dissipation in the under- 
lying fluid with the MD viscosity taking the place of a 
postulated friction coefficient. In addition there is a dis- 
sipative term parallel to U^; . The energy part of Eq. (^) 
is identical in form to the energy equation postulated 
by Avalos and Mackic and similar to the equation 
studied by Espahol ||T^, save for the fact that here the 
work done by the conservative force (pk + Pi)^ki ■ Ufc;/4 
is present. The principal difference is associated with 
the Voronoi lattice: Our DP's fill space and change their 
shape, a key feature that enables this new DPD to treat 
a multitude of length scales within a single simulation. 

In order to obtain F and q we invoke the Marko- 
vian approximation to write F = Wfejn IVfcjn -I- u}ki±Wki^, 
where Yki is decomposed into components parallel and 
perpendicular to e^;, and the PF's are defined as Gaus- 
sian random variables with the correlation function 

{Wkla{t)Wnml3{t')) = 5ap5{t -t'){5kn5lni+ ^kin^ln) whcrC 

a and (3 denotes either _L or ||. Newton's third law guar- 
antees that u>ki = —<^ik- We have the similar expression 
m = AfejW^fci where A^; = -Ajfc and Wki satisfies the 
above equation for the VF's without the Ja/j-factor. 

In Refs. 1^,^^ the magnitudes of Ffc; and qu 
are obtained on the basis of the Fokker-Planck equa- 
tion which derives from equations like Eqs. (^. The 
isothermal results, adapted to the present case, are 
the fluctuation dissipation relations a;^.;|| = 2ajf;^ ~ 

^r]kBT{lki/rki) for the force F^; and for the heat fluc- 
tuations A^; — 2kBT\{lki/rki) fl^. It is also pos- 
sible to show that detailed balance holds, and 
that the DP's obey the Gibbs distribution = 
Z-i(r,T/)exp(-/3Ejp2/2Mfe-HF(rfc))), where the 
potential V{Yk) is responsible for the pressure force in 
Eq. and T = l/{l3kB) is the temperature character- 
izing the MD system [|l2| . 

Considering now the application illustrated in Fig. ^ 
we need to define DP-coUoid forces. Taking the hydro- 
dynamic momentum flux tensor and Eq. (^) as a starting 
point we observe that the DP-coUoid interaction may be 
obtained in the same form as the DP-DP interaction of 
Eq. with the replacement of Iki Lki, where Lki 
is the length (area in 3D) of the arc segment where the 
DP meets the colloid (see Fig. |l|). The velocity gradi- 
ent is that between the DP and the colloid surface. The 
latter may be computed by linear interpolation using Ufc 
and the velocity of the colloid surface together with a no- 
slip boundary condition on this surface. In the momen- 
tum fluctuation-dissipation relation too the replacement 
Iki Lki must be made. In order to increase the spatial 
resolution where colloidal particles are close it is neces- 
sary to introduce a higher DP density there; this ensures 
that fluid lubrication effects are maintained. After these 
particles have moved it may be necessary to re-tile the 
DP system, as is done in finite element calculations. This 



is simply achieved by distributing the mass, momentum 
and energy of the old DP's on the new ones according to 
their area (or volume in 3D). When a new Voronoi cell is 
created on top of the old ones the new DP mass is sim- 
ply the sum of the old DP mass densities times the area 
of overlap between the new and old dissipative particles 
(and similarly for the DP momentum and energy). This 
technique may also be used to model polymer solutions 
if the polymer chains are formed of linked beads. 

The DPD which we have derived in the present work 
is similar to conventional DPD, without- |^ and with 
energy conservation [ p^ . But the forces conventionally 
used to define DPD have now been given a microscopic 
basis. More important, however, is the fact that our 
analysis permits the introduction of specific physical in- 
teractions at the mesoscopic level, together with a well- 
defined meaning for this mesoscale. Finally, we note the 
similarity of the present particulate description, which is 
based on a bottom-up approach, to existing continuum 
approaches, which start out from a macroscopic descrip- 
tion. Such top-down approaches include in particular 
smoothed particle hydrodynamics [p^ and finite-element 
simulations. In these descriptions too the computational 
method is based on tracing the motions of elements of the 
fluid on the basis of the forces acting between them . 
We stress, however, that while such top-down compu- 
tational strategies require as initial input a macroscopic 
phcnomcnological description, the present approach re- 
lies on a microscopic representation from the outset. 
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